home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / LAGUER.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  4KB  |  137 lines
  1. PROCEDURE laguer(a: glcarray; m: integer; VAR x: gl2array;
  2.        eps: real; polish: boolean);
  3. (* Programs using routine LAGUER must define in the main routine the types
  4. TYPE
  5.    glcarray = ARRAY [1..2*m+2] OF real;
  6.    gl2array = ARRAY [1..2] OF real; *)
  7. LABEL 99;
  8. CONST
  9.    epss=6.0e-8;
  10.    mxit=100;
  11. VAR
  12.    j,iter: integer;
  13.    err,dxold,cdx,abx,dum: real;
  14.    sq,h,gp,gm,g2,g: gl2array;
  15.    b,d,dx,f,x1,cdum: gl2array;
  16. PROCEDURE cdiv(a,b: gl2array; VAR c:gl2array);
  17. (* Complex division of a by b, answer in c *)
  18. VAR
  19.    r,den: real;
  20. BEGIN
  21.    IF (abs(b[1]) >= abs(b[2])) THEN BEGIN
  22.       r := b[2]/b[1];
  23.       den := b[1]+r*b[2];
  24.       c[1] := (a[1]+a[2]*r)/den;
  25.       c[2] := (a[2]-a[1]*r)/den
  26.    END ELSE BEGIN
  27.       r := b[1]/b[2];
  28.       den := b[2]+r*b[1];
  29.       c[1] := (a[1]*r+a[2])/den;
  30.       c[2] := (a[2]*r-a[1])/den
  31.    END
  32. END;
  33. FUNCTION cabs(a: gl2array): real;
  34. (* Complex absolute value of a *)
  35. VAR
  36.    x,y : real;
  37. BEGIN
  38.    x := abs(a[1]);
  39.    y := abs(a[2]);
  40.    IF (x = 0.0) THEN cabs := y
  41.    ELSE IF (y = 0.0) THEN cabs := x
  42.    ELSE IF (x > y) THEN cabs := x*sqrt(1.0+sqr(y/x))
  43.    ELSE cabs := y*sqrt(1.0+sqr(x/y))
  44. END;
  45. PROCEDURE csqrt(a: gl2array; VAR b: gl2array);
  46. (* Returns complex square root of a in b *)
  47. VAR
  48.    x,y,u,v,w,r : real;
  49. BEGIN
  50.    IF ((a[1] = 0.0) AND (a[2] = 0.0)) THEN BEGIN
  51.       u := 0.0;
  52.       v := 0.0
  53.    END ELSE BEGIN
  54.       x := abs(a[1]);
  55.       y := abs(a[2]);
  56.       IF (x >= y) THEN
  57.          w := sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+sqr(y/x))))
  58.       ELSE BEGIN
  59.          r := x/y;
  60.          w := sqrt(y)*sqrt(0.5*(r+sqrt(1.0+sqr(r))))
  61.       END;
  62.       IF (a[1] >= 0.0) THEN BEGIN
  63.          u := w;
  64.          v := a[2]/(2.0*u)
  65.       END ELSE BEGIN
  66.          IF (a[2] >= 0.0) THEN v := w
  67.          ELSE v := -w;
  68.          u := a[2]/(2.0*v)
  69.       END
  70.    END;
  71.    b[1] := u;
  72.    b[2] := v;
  73. END;   
  74. BEGIN
  75.    dxold := cabs(x);
  76.    FOR iter := 1 TO mxit DO BEGIN
  77.       b[1] := a[2*m+1];
  78.       b[2] := a[2*m+2];
  79.       err := cabs(b);
  80.       d[1] := 0.0;
  81.       d[2] := 0.0;
  82.       f[1] := 0.0;
  83.       f[2] := 0.0;
  84.       abx := cabs(x);
  85.       FOR j := m DOWNTO 1 DO BEGIN
  86.          dum := f[1];
  87.          f[1] := x[1]*f[1]-x[2]*f[2]+d[1];
  88.          f[2] := x[1]*f[2]+x[2]*dum+d[2];
  89.          dum := d[1];
  90.          d[1] := x[1]*d[1]-x[2]*d[2]+b[1];
  91.          d[2] := x[1]*d[2]+x[2]*dum+b[2];
  92.          dum := b[1];
  93.          b[1] := x[1]*b[1]-x[2]*b[2]+a[2*j-1];
  94.          b[2] := x[1]*b[2]+x[2]*dum+a[2*j];
  95.          err := cabs(b)+abx*err
  96.       END;
  97.       err := epss*err;
  98.       IF (cabs(b) <= err) THEN BEGIN
  99.          dx[1] := 0.0;
  100.          dx[2] := 0.0;
  101.          GOTO 99
  102.       END ELSE BEGIN
  103.          cdiv(d,b,g);
  104.          g2[1] := sqr(g[1])-sqr(g[2]);
  105.          g2[2] := 2.0*g[1]*g[2];
  106.          cdiv(f,b,cdum);
  107.          h[1] := g2[1]-2.0*cdum[1];
  108.          h[2] := g2[2]-2.0*cdum[2];
  109.          cdum[1] := (m-1)*(m*h[1]-g2[1]);
  110.          cdum[2] := (m-1)*(m*h[2]-g2[2]);
  111.          csqrt(cdum,sq);
  112.          gp[1] := g[1]+sq[1];
  113.          gp[2] := g[2]+sq[2];
  114.          gm[1] := g[1]-sq[1];
  115.          gm[2] := g[2]-sq[2];
  116.          IF(cabs(gp) < cabs(gm)) THEN BEGIN
  117.             gp[1] := gm[1];
  118.             gp[2] := gm[2]
  119.          END;
  120.          cdum[1] := m;
  121.          cdum[2] := 0.0;
  122.          cdiv(cdum,gp,dx);
  123.       END;
  124.       x1[1] := x[1]-dx[1];
  125.       x1[2] := x[2]-dx[2];
  126.       IF ((x[1] = x1[1]) AND (x[2] = x1[2])) THEN GOTO 99;
  127.       x[1] := x1[1];
  128.       x[2] := x1[2];
  129.       cdx := cabs(dx);
  130.       IF ((iter > 6) AND (cdx >= dxold)) THEN GOTO 99;
  131.       dxold := cdx;
  132.       IF (not polish) THEN
  133.          IF (cabs(dx) <= eps*cabs(x)) THEN GOTO 99
  134.    END;
  135.    writeln('pause in routine LAGUER - too many iterations'); readln;
  136. 99:   END;
  137.